package com.example.analysis.utils;

import java.util.Arrays;
import java.util.Random;
import static com.example.analysis.utils.ComCalcUtils.calcMean;
import static com.example.analysis.utils.ComCalcUtils.calcStd;
import static org.apache.commons.math3.special.Erf.erf;

public class NormalityTest {
    public static double alpha = 0.05;

    // ks检验。pValue > 0.05，呈现正态
    public static double[] ksTest(double[] data) {
        // 1.计算均值
        double mean = calcMean(data);
        // 2.计算标准差
        double std = calcStd(data, mean);
        System.out.println("均值：" + mean + "，标准差：" + std);
        // 3.计算D
        double d = calcStatistic(data, mean, std);
        // 4.计算P值
        double pValue = calcPValue(d, data, 1000);
        System.out.println("Lilliefors 检验：dValue = " + d + ", pValue = " + pValue);
        return new double[]{mean, std, d, pValue};
    }

    /**
     * 计算D值，通过比较经验分布函数（ECDF）和理论正态分布函数（CDF）的最大差距。
     * D越小，P越大，越能满足正态分布
     *
     * @param data  数据
     * @param mu    均值
     * @param sigma 标准差
     * @return D值
     */
    public static double calcStatistic(double[] data, double mu, double sigma) {
        if (sigma == 0)
            return 0;
        // 排序
        Arrays.sort(data);

        int n = data.length;
        double maxDiff = 0.0;

        for (int i = 0; i < n; i++) {
            double x = data[i];
            double erf = erf((x - mu) / (sigma * Math.sqrt(2)));
            double cdf = 0.5 * (1 + erf); // 理论 CDF
            double dPlus = (i + 1.0) / n - cdf;    // 正向偏离
            double dMinus = cdf - i / (double) n;  // 负向偏离
            maxDiff = Math.max(maxDiff, Math.max(dPlus, dMinus));
        }
        return maxDiff;
    }

    /**
     * 计算P值
     * 近似公式计算 p 值（Marsaglia 方法，适用于 n >= 30）
     *
     * @param d              D值
     * @param data           数据
     * @param numSimulations 模拟次数
     * @return 返回P值
     */
    public static double calcPValue(double d, double[] data, int numSimulations) {
        int n = data.length;
        if (n == 0)
            return 0;
        double mean = calcMean(data);
        double std = calcStd(data, mean);
        Random random = new Random();

        int count = 0;
        for (int i = 0; i < numSimulations; i++) {
            // 模拟数据
            double[] simulatedData = new double[n];
            for (int j = 0; j < n; j++) {
                // nextGaussian()用于生成服从标准正态分布的随机数。
                simulatedData[j] = random.nextGaussian() * std + mean;
            }
            double simulatedMean = calcMean(simulatedData);
            double simulatedStd = calcStd(simulatedData, simulatedMean);
            double simulatedD = calcStatistic(simulatedData, simulatedMean, simulatedStd);
            // 统计统计量大于样本统计量的个数
            if (simulatedD >= d) {
                count++;
            }
        }
        return (double) count / numSimulations;
    }
}
